In [1]:
%matplotlib inline
import pandas as pd
import zipfile
import shutil
import urllib2
from urllib2 import urlopen
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as dates
import matplotlib.ticker as tick
import scipy.stats as sp
import statsmodels.api as sm
from pandas.stats.api import ols
from datetime import datetime
from bs4 import BeautifulSoup
from pylab import rcParams
import platform
rcParams['figure.figsize'] = 15, 10
import re
import os
import sys
import glob
import urllib
import HTMLParser
from cStringIO import StringIO
import pyproj
from pyproj import Proj
import gzip
import ftplib
import calendar
import datetime
from datetime import date
import pymodis
In [2]:
import arcpy
arcpy.CheckOutExtension("spatial")
from arcpy import env
from arcpy.sa import *
In [3]:
print("Operating System " + platform.system() + " " + platform.release())
print("Python Version " + str(sys.version))
print("Pandas Version " + str(pd.__version__))
print("Numpy Version " + str(np.__version__))
In [4]:
import wellapplication as wa
import UBM
In [5]:
Zonal_HUCS = "H:/GIS/BCM/ZonalExtra.gdb/Zonal_HUC"
Zone_field = "HUC_12"
z_Name = "H:/GIS/BCM/ZonalExtra.gdb"
arcpy.env.overwriteOutput = True
In [ ]:
indata = "H:/GIS/BCM/AvailableWater.gdb"
UBM.zone_gdb(indata, z_Name, Zonal_HUCS, Zone_field)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
def replace_hdr_file(hdrfile):
"""
Replace the .hdr file for a .bil raster with the correct data for Arc processing
Required: hdrfile -- filepath for .hdr file to replace/create
Output: None
"""
# hdr file replacment string
HDRFILE_STRING = "byteorder M\n\
layout bil\n\
nbands 1\n\
nbits 16\n\
ncols 6935\n\
nrows 3351\n\
ulxmap -124.729583333331703\n\
ulymap 52.871249516804028\n\
xdim 0.00833333333\n\
ydim 0.00833333333\n"
with open(hdrfile, 'w') as o:
o.write(HDRFILE_STRING)
In [ ]:
path= "H:/GIS/MODIS2/MODIS.gdb/"
arcpy.env.workspace = path
for rast in arcpy.ListRasters()[:3]:
dsc = arcpy.Describe(rast)
print(dsc.baseName)
In [ ]:
for hdrfile in glob.glob("H:/GIS/SNODAS/SNODASUNZIPPED/SNODASUNGZ/*.Hdr"):
replace_hdr_file(hdrfile)
In [ ]:
sr = arcpy.SpatialReference(4326) #Spatial Projection WGS84
# Set output coordinate system
outCS = arcpy.SpatialReference('NAD 1983 UTM Zone 12N') #NAD83 Zone 12 Code is 26912
# Set local variables
inMaskData = "H:/GIS/SNODAS/SNODAS.gdb/UT_HUC_area"
# from Table 4 of http://nsidc.org/data/docs/noaa/g02158_snodas_snow_cover_model/
prodcodelong = {1025: 'Precipitation', 1034: 'Snow water equivalent', 1036: 'Snow depth',
1038: 'Snow pack average temperature', 1039: 'Blowing snow sublimation',
1044: 'Snow melt', 1050: 'Snow pack sublimation'}
prodcode = {1025: 'PREC', 1034: 'SWEQ', 1036: 'SNOD', 1038: 'SPAT', 1039: 'BSSB', 1044: 'SNML', 1050: 'SPSB'}
# iterate through data files creating rasters and defining projections
for dtfile in glob.glob("H:/GIS/SNODAS/SNODASUNZIPPED/SNODASUNGZ/*.dat"):
indata = arcpy.Raster(dtfile)
arcpy.DefineProjection_management(indata, sr)
typAbb = prodcode[int(os.path.basename(dtfile)[8:12])]
#indata.save("H:/GIS/SNODAS/SNODAS.gdb/" + typAbb + dtfile[-20:-11]) #YYYYMMDD
# Execute ExtractByMask to clip snodas data to Utah watersheds
outExtractByMask = arcpy.sa.ExtractByMask(indata, inMaskData)
# Determine the new output feature class path and name = productcode + YYYYMMDD and save to a geodatabase
outrast = "H:/GIS/SNODAS/SNODAS.gdb/" + typAbb + dtfile[-20:-11] #os.path.join(outWorkspace, rast)
# Print name to veryfy save
print(typAbb + dtfile[-20:-11])
# Project Raster to UTM Zone 12
arcpy.ProjectRaster_management(outExtractByMask, outrast, outCS, 'BILINEAR', '1000',\
'WGS_1984_(ITRF00)_To_NAD_1983', '#', '#')
From http://nsidc.org/data/docs/noaa/g02158_snodas_snow_cover_model/, the file abbreviations are as follows:
In [ ]:
prodcode = {'us_ssmv11038wS__A':'SPAT', 'us_ssmv11044bS__T':'SNML', 'us_ssmv11050lL00T':'SPSB',
'us_ssmv11034tS__T':'SWEQ', 'us_ssmv01025SlL00':'RAIN', 'us_ssmv01025SlL01':'SNOW',
'us_ssmv11036tS__T':'SNOD', 'us_ssmv11039lL00T':'BSSB'}
path = "H:/GIS/SNODAS/geotifSNODAS/SNDS/"
for filename in os.listdir(path):
if filename.startswith("us_ssmv"):
code = prodcode[filename[0:17]]
yrsrt = filename.find('TNATS') + 5
yr = filename[yrsrt:yrsrt+4]
mo = filename[yrsrt+4:yrsrt+6]
dy = filename[yrsrt+6:yrsrt+8]
os.rename(os.path.join(path, filename), os.path.join(path,code+yr+mo+dy+filename[-4:]))
In [ ]:
def mergeRasts(path, data_type = 'AET', monthRange = [1,4], yearRange = [2000,2001]):
arcpy.env.workspace = path
print(arcpy.ListRasters())
for y in range(yearRange[0],yearRange[-1]+1): #set years converted here
for m in range(monthRange[0],monthRange[-1]+1): #set months converted here
nm = data_type + str(y) + str(m).zfill(2)
rlist=[]
for rast in arcpy.ListRasters(nm+'*'):
rlist.append(rast)
print(rlist)
In [ ]:
path="H:/GIS/MODIS2/MODIS.gdb"
mergeRasts(path)
In [ ]:
path = "H:/GIS/SNODAS/geotifSNODAS/SNDS/"
g = {}
from arcpy.sa import *
arcpy.env.workspace = path
arcpy.env.overwriteOutput = True
code = 'PREC'
for y in range(2003,2016):
for m in range(1,13):
g[code+str(y)+str(m).zfill(2)] = []
for name in sorted(glob.glob(path+code+'*.tif')):
rast = os.path.basename(name)
if rast[0:4] == code and int(rast[4:8]) == y and int(rast[8:10]) == m:
g[code+str(y)+str(m).zfill(2)].append(rast)
else:
pass
if len(g[code+str(y)+str(m).zfill(2)])>0:
print(g[code+str(y)+str(m).zfill(2)])
calc = CellStatistics(g[code+str(y)+str(m).zfill(2)], statistics_type = "SUM", ignore_nodata="DATA")
calc.save("H:/GIS/SNODAS/SNODAS.gdb/"+rast[0:4]+str(y).zfill(2)+str(m).zfill(2)+"SUM")
In [ ]:
In [ ]:
sr = arcpy.SpatialReference(4326) #Spatial Projection WGS84
# Set output coordinate system
outCS = arcpy.SpatialReference('NAD 1983 UTM Zone 12N') #NAD83 Zone 12 Code is 26912
# Set local variables
inMaskData = "H:/GIS/SNODAS/SNODAS.gdb/UT_HUC_area"
# from Table 4 of http://nsidc.org/data/docs/noaa/g02158_snodas_snow_cover_model/
prodcodelong = {1025: 'Precipitation', 1034: 'Snow water equivalent', 1036: 'Snow depth',
1038: 'Snow pack average temperature', 1039: 'Blowing snow sublimation',
1044: 'Snow melt', 1050: 'Snow pack sublimation'}
prodcode = {1025: 'PREC', 1034: 'SWEQ', 1036: 'SNOD', 1038: 'SPAT', 1039: 'BSSB', 1044: 'SNML', 1050: 'SPSB'}
# iterate through data files creating rasters and defining projections
for dtfile in glob.glob("H:/GIS/SNODAS/SNODASUNZIPPED/SNODASUNGZ/*.dat"):
indata = arcpy.Raster(dtfile)
arcpy.DefineProjection_management(indata, sr)
typAbb = prodcode[int(os.path.basename(dtfile)[8:12])]
# Determine the new output feature class path and name = productcode + YYYYMMDD and save to a geodatabase
outrast = "H:/GIS/SNODAS/SNODAS.gdb/" + typAbb + dtfile[-20:-11] #os.path.join(outWorkspace, rast)
# Print name to veryfy save
print(typAbb + dtfile[-20:-11])
# Project Raster to UTM Zone 12
arcpy.ProjectRaster_management(outExtractByMask, outrast, outCS, 'BILINEAR', '1000',\
'WGS_1984_(ITRF00)_To_NAD_1983', '#', '#')
In [ ]:
arcpy.env.overwriteOutput = True
path="H:/GIS/SNODAS/SNODASproj.gdb/"
arcpy.env.workspace = path
for name in arcpy.ListRasters():
In [ ]:
monthRange = [1,12]
yearRange = [2003,2016]
g = {}
path="H:/GIS/SNODAS/SNODASproj.gdb/"
arcpy.env.workspace = path
arcpy.env.overwriteOutput = True
for y in range(yearRange[0],yearRange[1]+1): #set years converted here
for m in range(monthRange[0],monthRange[1]+1): #set months converted here
my = str(y)+str(m).zfill(2)
newdn = 'TPPT' + my
try:
calc = Plus('SNOW'+ my +'SUM', 'RAIN'+ my +'SUM')
calc.save(newdn+'SUM')
print(newdn)
except(RuntimeError):
pass
In [ ]:
from arcpy.sa import *
monthRange = [1,12]
yearRange = [2005,2005]
path = "H:/GIS/Calc.gdb/"
path1 = "H:/GIS/SNODAS/SNODASproj.gdb/"
path2 = "H:/GIS/MODIS/MODOUT.gdb/"
arcpy.env.workspace = path
arcpy.env.overwriteOutput = True
area = 'H:/GIS/NHD_UT_Proj.gdb/UT_HUC_area'
arcpy.env.mask = area
for y in range(yearRange[0],yearRange[1]+1): #set years converted here
for m in range(monthRange[0],monthRange[1]+1): #set months converted here
my = str(y)+str(m).zfill(2)
newdn = 'AVLW' + my
rain = (path1 + 'RAIN'+ my +'SUM')
snowMelt = (path1 + 'SNML'+ my +'SUM')
actEvap = (path2 + 'AET'+ my)
avail = (Con(IsNull(rain),0, rain)) + (Con(IsNull(snowMelt),0, snowMelt)) - (Con(IsNull(actEvap),0, actEvap))
avail = (available < 0, 0, available)
available.save(newdn)
print(newdn)
In [ ]: